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SUMMARY 


Under  the  Office  of  Naval  Research(ONR)  Grant  No.N00014-84-K-0610 ,  the 
University  of  the  District  of  Coliiinbia(UDC)  has  undertaken  a  basic  research  on 
dynamics  of  marine  atmospheres  in  tropical  cyclones  with  an  objective  of 
improving  understanding  of  turbulence  in  cyclone  vortices  and  efficiency  of 
computer  simulation  for  transfer  of  heat,  moisture  and  momentum  in  marine 
planetary  boundary  layers  of  tropical  cyclones.  For  this  objective,  UDC  has 
performed  the  following  four  phases  of  effort; 

I.  Develop  a  mathematical  model,  using  the  second 

order  turbulent  closure  eguations,  for  simulating 
planetary  boundary  layers  of  tropical  cyclones 
over  the  sea  surface. 


II.  Calculate  velocity,  enthalpy  and  moisture  dis¬ 
tributions  in  planetary  boundary  layers  of  tropical 
cyclones  over  the  sea  surface. 

III.  Develop  and  instrument  a  laboratory  vortex 
chamber  for  measuring  the  air/water  interface 
turbulence  and  heat,  moisture  and  momentum 
transfer . 

IV.  Develop  a  graphics-oriented  interactive 
finite-element  computer  model  to  simulate 
transfer  of  heat,  moisture  and  momentum 
in  marine  planetary  boundary  layers  of 
tropical  cyclones. 

In  this  report,  results  of  these  four  phases  of  effort  are  presented;  and 


conclusions  of  this  research  program  and  recommendations  for  further  research 
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I .  INTRODUCTION 

The  structure,  evolution  and  motion  of  tropical  cyclones  are  controlled 
by  complex  interaction  of  microscale,  convective,  mesoscale  and  synoptic 
processes.  Convergence  of  moisture  in  the  planetary  boundary  layer  is  a 
primary  mechanism  for  organization  of  the  convection  in  the  inner  core  of  an 
intense  tropical  cyclone.  The  convection,  in  turn,  provides  the  diabatic 
heating  necessary  to  sustain  intensity  of  the  tropical  cvclone.  Recent 
advances  in  turbulent  modeling  and  computational  techniques  have  made  it 
possible  to  simulate  simultaneous  effects  of  numerous  parameters  on  dynamics 
of  tropical  cyclones. 

Under  the  Office  of  Naval  Research(ONR)  Grant  No.  N00014-84-K-0610 ,  the 
University  of  the  District  of  Columbia(UDC)  has  undertaken  a  basic  research  on 
dynairdcs  of  marine  atmospheres  in  tropical  cyclones  with  an  objective  of 
improving  understanding  of  turbulence  in  cyclone  vortices  and  -efficiency  of 
computer  simulation  for  transfer  of  heat,  moisture  and  momentum  in  marine 
planetary  boundary  layers  of  tropical  cyclones.  For  this  objective,  UDC  has 
performed  the  following  four  phases  of  effort: 

I.  Develop  a  mathematical  model,  using  the  second 
order  turbulent  closure  equation,  for  simulating 
planetary  boundary  layers  of  tropical  cyclones 
over  the  sea  surface. 

II.  Calculate  velocity,  enthalpy  and  moisture  dis¬ 
tributions  in  planetary  boundary  layers  of  tropical 
cyclones  over  the  sea  surface. 

III.  Develop  and  instrument  a  laboratory  vortex 
chamber  for  measuring  the  air/water  interface 
turbulence  and  heat,  moisture  and  momentum 
transfer. 

IV.  Develop  a  graphics-oriented  interactive 
finite-element  computer  model  to  simulate 
transfer  of  heat,  moisture  and  momentum  in 
marine  planetary  boundary  layers  of  tropical 
cyclones . 
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In  this  report,  a  mathematical  model,  using  the  second  order  turbulent 
closure  equations,  for  simulating  planetary  boundary  layers  of  tropical 
cyclones  over  the  sea  surface  is  described  in  Section  2.  Results  of  calcu¬ 
lated  velocity,  enthalpy  and  moisture  distributions  in  planetary  boundary 
layers  of  tropical  cyclones  over  the  sea  surface  are  presented  in  Section  3. 
To  increase  confidence  in  the  theory,  calculated  results  are  compared  with 
laboratory  test  data  in  section  4.  Development  of  a  graphics-oriented  inter¬ 
active  finite-element  computer  program  for  simulating  the  vortex  flow  consti¬ 
tutes  a  major  effort  in  a  latter  part  of  this  research  program;  and  it  will  be 
detailed  in  Section  5  of  this  report.  Conclusions  of  this  research  effort  and 
recommendations  for  further  research  on  dynamics  of  marine  atmospheres  in 
tropical  cyclones  are  given  in  Section  6. 

2.  A  MATHEMATICAL  MC»EL  FOR  SIMULATING 

VORTEX  PLANETARY  BOUNDARY  LAYERS 

Figure  1  shows  schematically  a  vortex  marine  planetary  boundary  layer. 

The  situation  presents  many  interesting  features  which  include  the  boundary 

layer  heat,  moisture  and  momentum  transfer;  turbulent  energy  production, 

dissipation  and  diffusion;  strong  influence  of  the  gravitational,  Coriolis  and 

centrifugal  forces;  and  influences  of  the  air/sea  interface  roughness  and 

liquid  droplets  entrainment.  Under  Phase  I  effort  of  this  research  program,  a 
\> 

mathematical  model,  using  the  second  order  turbulent  closure  equations,  for 
simulating  planetary  boundary  layers  of  tropical  cyclones  over  the  sea  surface 
has  been  developed. 

Previous  discussions  of  the  turbulent  transport  theory  (e.g.,  see 
Harlow  and  Hirt,  1969;  Mellor  and  Yamada,  1974;  and  Leweller,  1981)  yielded 
equations  in  rectangular  coodinates  for  the  turbulent  stresses  components  and 
flux  moments.  For  simulating  the  vortex  marine  planetary  boundary  layers,  it 
is  more  convenient  to  use  the  cylindrical  coodinates.  Under  Phase  I  effort  of 
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this  program,  a  second-order  turbulent  transport  theory  which  uses  the 
cylindrical  coordinates  and  accounts  for  special  features  which  are  described 
in  the  preceding  paragraph  is  derived.  Details  of  derivation  of  the  theory 
are  described  in  Chi  (1986);  but  a  brief  summary  of  the  theory  and  its 
resultant  equations  are  presented  below  in  this  section. 

Starting  from  instantaneous  conservation  equations,  mean  equations  of 
motion,  and  balances  of  energy  and  moisture  for  turbulent  vortex  flows  were 
first  derived.  Transport  equations  for  turbulent  stress  components,  and 
enthalpy  and  moisture  flux  moments  were  then  derived.  Closure  assumptions 
formulated  by  Chou  (1945),  Launder,  Reese  and  Rodi  (1975),  and  Lumley  and 
Newman  (1977)  were  finally  used  to  complete  the  closure.  In  conclusion,  a 
mathematical  model  for  simulating  vortex  planetary  boundary  layers  has  been 
developed . 

The  coordinates  are  chosen  so  that  z  axis  is  vertical,  and  r  axis  is  in 
the  radial  direction  and  4  is  the  azimuth  angle.  Equations  governing  mean 
motion  and  balances  of  mean  enthalpy  and  moisture  concentration  can  be  written 
as ; 
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In  these  equations,  u,  v  and  w  are  the  mean  velocity  components  in  the  r, 

<}  and  z  directions;  t  is  the  time,  p  the  mean  pressure,  h  the  mean  moist 

specific  enthalpy  which  is  defined  as  (s+Lq^),  and  q  the  mean  total  moisture 

mixing  ratio  which  is  defined  as  (q^+qj^).  Here,  s  is  the  mean  dry  specific 

enthalpy  which  is  defined  as  Cp(T-(Te-g2/Cp) ] ,  L  the  water  latent  heat  of 

vaporization,  the  constant  pressure  specific  heat,  T  the  temperature,  Tc 

the  standard  temperature,  the  water  vapor  mixing  ratio,  q^^^  the  liquid  water 

mixing  ratio,  g  the  gravitational  acceleration,  f  the  Coriolis  parameter, 

the  virtual  dry  specific  enthalpy  which  is  defined  as  s( 1+1 .609Qq^.  -q),  Cc  the 

standard  density,  6  the  thermal  expansion  coefficient  which  is  defined  as 

(l/(CpT. )],  the  turbulent  kinematic  viscosity,  and  the  turbulent  Prandtl 

number,  and  c,  the  turbulent  Schmidt  number, 
n 

In  equations  1  through  6,  the  turbulent  viscosity  value  can  be 

calculated  from  the  turbulent  kinetic  energy  k  and  dissipation  e  values  as 
follows ; 
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In  quations  (8)  and  (9),  the  virtual  dry  specific  enthalpy  value  value  is 
related  to  the  mean  moist  specific  enthalpy  h  and  mean  total  moisture  mixing 
ratio  q  values  by  the  equations: 
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where  y  and  t  are  equal  to  (L/C  )(3q/3T)  and  (C  T/L),  respectively. 

P  P 

Equations  (1)  through  (11)  with  appropriate  initial  and  boundary 
conditions  now  constitute  a  ccHtiplete  equation  set  which  may  be  solved  to 
simulate  marine  planetar^  boundary  layers  of  tropical  cyclones. 

3.  HEAT,  MOISTURE  AND  MOMENTUM  DISTRIBUTIONS 
IN  VORTEX  BOUNDARY  LAYERS 

The  turbulent  models  discussed  in  Section  2  above  can  be  used  to  simulate 
the  marine  planetary  boundary  layer  of  tropical  cyclones.  As  an  example,  a 
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simple  case  which  will  first  be  considered  is  that  the  flow  is  quasi-steady 
and  axisymmetrical  and  that  the  boundary  layer  thickness  is  much  less  than  the 
vortex  radial  length  (i.e.,  for  a  region  of  the  vortex  boundary  layer  far  away 
from  the  vortex  center).  Governing  equations  for  this  problem  can  be  derived 
from  equations  which  are  discussed  in  Section  2  above  and  they  may  be  written 
as : 
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where 


r 

V  =  C  — 

t  V  E 


Ih  =  V  c  -  4  ^ ^  -^]  (For  Clear  Region) 


3z  t  6  c.  Sz 

n 


0  dZ 

q 


—  =  V  c  ^  1-609 iil  +  IL  13.]  (for  Cloudy  Region) 
3z  t  e  (1  +  y)c.  3z  0  3z' 


Cv  =  0.09exp[-2.5/(l+0,02Re^) 


Cj,  -  2[1  -  0 . 3exp(-Re^ )  ] 


Re^  =  — 
t  e 


(c,  ,  0  ,  0  ,  0  ,  0  )  =  (0.9,  0.9,  0.9,  1.0,  1.3) 

w  ^  6  K  C 


(19) 

(20) 

(21) 

(22) 

(23) 

(24) 

(25) 


A  finite  difference  numerical  procedure  has  been  developed  and  a  computer 
program  has  been  written  for  solving  the  above  set  of  equations  with 
appropriate  vortex  boundary  conditions  to  simulate  the  transfer  of  heat, 
moisture  and  momentum  in  turbulent  vortex  flows  over  the  water  surface. 
Details  of  this  development  has  been  published  by  Chi  (1987),  and  a  reprint  of 
the  paper  is  appended  to  this  report  (see  Appendix  A). 

As  an  example,  the  coir^juter  results  for  a  hypothetic  tropical  cyclone 
under  conditions  shown  in  table  1  are  presented  and  discussed  below: 

TABLE  1 

Maximum  tangential  Wind  v+  at  1  km  Alt: tude  52  m/s 
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Radius  at  Which  Max.  v+  Occurs 


40  km 


Temperature  at  1  km  Altitude  273  K 

Coriolis  Force  Coefficient  1x10  ^  s  ^ 

Temperature  at  Sea  Surface  288  K 

Fig-ure  2  shows  distributions  of  radial  and  tangential  components  (u_^  and 
)  of  the  cyclone  mean  velocity  at  the  1-km  altitude.  Radial  and  tangential 
components  (u  and  v)  of  the  boundary- layer  velocity  are  plotted  versus  height 
z  at  several  radii  r  in  figure  3  and  4,  respectively.  Values  for  the 
vertical-velocity  component  at  the  1-km  altitude  are  plotted  in  figure  5. 
The  boundary- layer  specific  enthaly  (h-h^)  and  moisture  ratio  (q-q^ )  values 
are  plotted  versus  height  z  at  several  radii  in  figures  6  and  7,  respectively. 
Figure  8  show  a  distribution  of  the  boundary- layer  turbulent  kinetic  energy. 

Several  interesting  marine  vortex  boundary  layer  characteristics  can  be 
observed  from  these  figures.  For  the  region  far  above  the  ground  (i.e.,  at 
high  altidude),  it  can  be  seen  in  figure  4  that  tangential  velocity  dominates. 
The  balance  of  force  for  this  region  is,  therefore,  a  balance  of  pressure 
gradient  and  combination  of  the  centrifugal  and  Coriolis  forces.  Near  the  sea 
surface,  retardation  of  the  tangential  velocity  can  be  seen  in  figure  4.  This 
retardation  is  accompanied  by  a  reduction  in  centrifugal  and  Coriolis  forces. 
The  balance  of  pressure,  centrifugal  and  Coriolis  forces  is  thereby  destroyed. 
The  flow  in  this  region  is  then  characterized  by  entrainment  of  the  flow  into 
the  boundary  layer  as  indicated  by  the  large  induced  radial  velocity  shown  in 
figure  3  and  the  induced  downward  flow  (i.e.,  negative  w_^  values)  shown  in 
figure  5.  Figures  6  and  7  show  large  enthalpy  and  moisture  values  at  vicinity 
of  the  sea  surface  and  their  decreased  values  at  large  altitude.  Advection 
diffusion  of  the  turbulent  energy  in  the  boundary  layer  can  be  observed  from 
distribution  of  k  values  shown  in  figure  8. 

Detailed  measurement  of  the  boundary- layer  values  in  tropical  cyclone  is 
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*  tn\  ^  * 


FIGURE  M  COnPUTER  SIHULATED  TANGENTIAL  VELOCITY  DISTRIBUTION 
(V)  IN  THE  BOUNDARY  LAYER  OF  A  TROPICAL  CYCLONE 
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FIGURE  t  COnPUTER  SinULATED  HOIST  STATIC  ENTHALPY  DISTRIBUTION 

IN  THE  BOUNDARY  LAYER  OF  A  TROPICAL  CYCLONE 
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FIGURE  7  COMPUTER  SIMULATED  TOTAL  MOISTURE  RATIO  DISTRIBUTION 
(q-q^)  IN  THE  BOUNDARY  LAYER  OF  A  TROPICAL  CYCLONE 
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difficulty  if  not  impossible.  A  laboratory  vortex  chamber  with  air/water 
interface  has  been  built  at  UDC.  Some  UDC  test  data  are  presented  and  dis¬ 
cussed  in  the  next  section  of  the  report. 


4.  COKPARISCW  OF  CALCOLATED  RESULTS  WITH 

LABORATORY  TEST  DATA 

In  order  to  establish  confidence  in  the  ccwiputer  simulated  results, 
measurements  of  the  turbulent  vortex  flow  have  been  made  on  a  UDC  vortex 
chamber.  Figure  9  shows  a  schematic  diagram  of  the  UDC  vortex  chamber.  The 
chamber  overall  dimension  is  55-cm  diameter  by  90-cm  high.  It  generates  an 
air  vortex  of  30.48-cm  diameter  by  14-cm  high.  The  vortex  is  generated  by 
forcing  air  through  24  evenly  spaced  vanes  of  15-cm  high  placed  on  a  30.48-cm 
diameter  pitch  circle.  The  swirling  air  discharges  from  the  chamber  through  a 
3.81-cm  diameter  hole  located  at  center  of  the  chamber's  top  disk.  Bottom  of 
the  swirling  air  is  in  contact  with  a  pool  of  water  which  can  be  maintained  at 
a  constant  depth  and  at  desired  temperatures.  The  chamber  is  instrumented 
with  both  the  TSI  hot-film  anemometer  and  laser  doppler  velocimeter.  An  IBM 
PC  has  been  used  for  online  data  acquisition  and  processing. 

Figure  10  shows  an  example  of  the  measured  radial  u_^  and  tantengial  v^ 
components  of  air  velocity  in  the  main  vortex  flow  outside  the  boundary  layer. 
Also  can  be  seen  in  this  figure  is  that  u^  and  v^  values  for  this  vortex  flow 
can  be  correlated  accurately  by  the  equations: 


u+ 


{u'+.)  r  2.506v  j 

€0  •  tQ 


1 .4v  r  ^2 

-  - ^{1  -  exp(-1.253(-^)  J} 

m 


(26) 


(27) 
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where  suffix  +  represents  values  at  vortex  main  flow,  suffix  represents 
values  at  large  radius,  and  suffix  m  represents  values  at  a  position  in  the 
main  flow  where  the  vortex  tangential  velocity  is  maximum. 

Measurements  in  the  boundary  layer  have  so  far  been  made  only  with  the 
hot-film  probe  for  the  air  flow  over  solid  end  walls.  Results  of  the  measured 
radial  u  and  tangential  v  values  inside  the  vortex  boundary  layers  on  the 
smooth  and  rough  end  walls  are  shown  in  figures  11  and  12,  respectively.  For 
comparison,  computer  simulated  results  under  same  conditions  as  experiments 
are  also  plotted  in  these  figures.  Excellent  agreement  between  theory  and 
experiments  can  be  seen  in  these  figures. 

5.  A  GRAPHICS-ORIENTED  INTERACTIVE  SIMULATION  MODEL 

FOR  TRANSFER  OF  HEAT,  MOISTURE  AND  MOMENTUM 

The  last  phase  of  effort  under  this  program,  was  spent  in  developing  a 
graphics  interactive  turbulent  simulation  computer  program  based  upon  the 
Galerkin's  finite  element  method.  A  peculiar  feature  of  the  finite-element 
method  is  the  division  of  domain  into  a  finite  number  of  elements.  A  set  of 
basis  functions  is  locally  defined  over  each  element  so  that  one  only  needs  to 
consider  the  governing  equations  for  individual  elements.  After  the  element 
equations  are  obtained,  a  simple  assembly  of  these  equations  over  all  the 
elements  yields  the  global  system  of  equations.  Since  the  shape  and  size  of 
individual  elements  are  very  flexible,  local  refinement  of  resolution  can  .be 
easily  obtained  for  regions  of  strong  gradients.  Irregularly  shaped  domains 
can  also  be  easily  divided  into  elements  without  coding  difficulties.  These 
flexibilities  of  the  finite-element  method  make  it  especially  useful  for 
modeling  heat,  moisture  and  momentum  transport  in  the  marine  boundary  layers 
for  which  mesh  nesting  and  complex  interactions  are  important.  In  addition, 
ability  of  discretization  allows  easy  interfacing  of  the  finite  element 
computer  program  with  the  interactive  graphics  software.  The  author  and  his 


21 


FIGURE  I?  COnPARISION  OF  TEST  DATA  WITH  COHPUTER  SinULATF.^ 
BOUNDARY  TANGENTIAL  VELOlITY  (v)  VALUES 


co-workers  have  been  developing  under  this  Grant  a  finite  element  analysis 
computer  program  for  simulating  three  dimensional  tranportation  of  momentum, 
moisture  and  enthalpy,  and  have  so  far  obtained  solutions  for  the  one-  and 
two-dimensional  problems.  In  addition,  the  UDC's  finite-element  computer 
programs  have  been  interfaced  with  the  PDA's  PATRAN  graphics-oriented  software 
package.  This  interface  has  allowed  PATRAN  to  be  used  for  model  creation  and 
results  presentation  under  an  interactive  graphics-oriented  environment  and 
the  UDC's  finite  element  computer  program  to  be  used  for  analysis. 

For  convenience  of  documentation,  development  todate  of  the  graphics- 
oriented  interactive  finite  element  computer  program  at  UDC  is  presented 
below  in  several  sub-sections  as  follows: 

.  Equations  Under  Consideration 

Finite  Element  Formulation 

.  Solution  Procedures  and  Computer  Programming 

.  Interfacing  of  the  UDC  Finite-Element  Programs  with 

the  Graphics-Oriented  PATRAN  Modeling  Pre-  and  Post- 
Processor 

.  Results  and  Discussion. 

5.1  Equations  Under  Consideration 

The  equations  under  consideration  are  conservation  equations  derived  in 
Section  2  of  the  report.  Limiting  discussions  in  this  section  to  problems  of 
momentum  transfer  in  two  dimensions  eind  assuming  constant  turbulent  viscosity, 
the  following  governing  equations  in  nondimens ional  tensor  form  may  be  ob¬ 
tained; 


3U. 

M(U.)  --~  +  -“(U.U.  ♦-—Pi.. 

1  3t  3Xj  1  j  p,  ij 


3U. 

57  sT-  ''Vh  * 

J 


=  0 


(28) 
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B(U.)  =  a  u.  a  +  b. 

1  1  1  Re  ax.  j  1 

J 


where  suffices  i  and  j  follow  the  usual  tensor  rules,  and  their  values  are 
either  one  or  two  for  reference  to  nondimensional  quantities  in  the  x  or  y 
direction.  The  nondiinensional  quantities  in  these  equations  in  terms  of 
reference  length  and  velocity  v  are  defined  as  follows: 


Nondimensional  Corrilis  Component,  F= 


Nondimensional  Pressure,  P  =  — 

C  V 

0  r 


Nondimensional  Reynolds  Number,  Re  = 


Nondimensional  Velocity  in  x-Direction,  U  =  U,  =  — 

^  1  v 


Nondimensional  Velocity  in  y-Direction,  V  = 

r 


Nondimensional  Distance  in  x-Director,  X  =  =  — 

r 


Nondimensional  Distance  in  y-Direction,  Y  =  X^  = 


Nondimensional  Time,  t  ■ 

r 
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5.2  Finite  Element  Formulation 

Construction  of  a  Galerkin  finite  element  algorithm  is  now  to  be  made  to 
solve  for  the  velocity  components  U1(=U)  and  U2(=V)  from  equation  (28).  The 
flew  domain  to  be  analyzed  is  divided  into  small  subdomains  call  finite  ele¬ 
ments.  Governing  equations  are  then  transformed  into  finite  element  equa¬ 
tions.  Finally,  the  finite  element  equations  are  assembled  into  a  global 
system  of  equations  which  can  be  solved  to  yield  simulation  results  for  the 
problem. 

In  the  present  analysis;  rectangular  finite  elements  are  used  to  discreti¬ 
zing  the  flow  domain.  The  rectangular  elements  have  four  nodes.  The  inter¬ 
polating  equations  for  velocity  in  an  element  shown  in  figure  13  can  be 
written  as: 


U^(r,  i.Tij''  =  N(r  j  ,  n  2 )  {UI(t  ) ) 

where 

^(1  -  n )) ( 1  -  n 2)'', 

I 

(1  +  rij)  (  1  -  n;)  i 

N(^i  ,^2 )  =  ^ 

(1  +  rij)(  1  +  ri2)  . 

U  1  -  n  i)  (  1  ri;)/ 

and  ni=Xl/a=X/a  and  n2=x2/b=Y/b. 


(39) 


(AO) 
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FI6URE  13  SKETCH  OF  A  RECTANGULAR  FINITE-ELEMENT 
DOMAIN 
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Equation  (39)  expresses  U  and  V  values  in  terms  of  natural  coordinates 
and  nodal  values  UI  which  is  equal  to  U1  or  U2.  Construction  of  finite 
element  solution  algorithm  for  equation  (28)  subjects  to  several  constraints. 
It  is  firstly  required  to  satisfy  the  boundary  condition  equation  (30), 
secondly  to  equate  gradient  of  the  convection  part  of  equation  (28)  to  zero 
for  controlling  nonlinear ly  induced  instability  and  dispersion  error, i.e.. 


^  W  V  .  W  w  . 

=  0 

3X .  eT  j  X . 


(41) 


and  finally  to  satify  the  penalty  function  which  equates  pressure  to  the 
continuiy  function  define  by  equation  (29),  i.e.: 


3U  . 

r(p)  =  p  +  X— (42) 

OA  . 

J 


Application  of  the  Galerkin  weighted  residue  method  to  governing  equation 
(28)  for  individual  finite  elements  with  constraints  of  equations  (30),  (41) 
and  (42),  and  substitution  of  discretization  equation  (39)  can  now  be  made  to 
yield  the  element  algorithm  statement: 


{N}M(U.)dr/  -  e.;„  {N}-|-M(U.  )d.Q  + /_  {N}B(U.)dr 

W  X  1i(dA*X  1  1 

e  e  j  e 

+  Tl-{N}n(p)da  =  0  (43) 

W  0  A  . 

e  1 

The  resultant  element  contributions  can  be  assenioled  into  a  global 
system  of  ordinary  differential  equations  with  respect  to  the  nondimens ional 
time  T  as  follows: 


[CI]{UI}'  +  [KI]{UI}  +  [KIJKUJ)  +  (SI)  =  0 


(44) 
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Here,  [Cl]  is  the  capacity  matrix,  [KI]  represents  the  convection  operator  in 
equation  (28),  [KIJ]  provides  the  complete  cross-coupling  between  the  nodel 
{UI}  and  {UJ}  values,  and  {SI}  is  any  source  term. 

5.3  Solution  Procedure  and  Ccmputer  Program  UDCFLOW 

An  implicit  finite-difference  integration  procedure  is  used  to  evaluate 
numerically  the  ordinary  differential  equation  set  (44).  With  ©,  which  is 
within  zero  and  one,  as  an  implicit  parameter,  the  differential  equation  set 
(44)  can  be  approximated  by  the  algebraic  equation  set: 


FI  =  [Cl]  UI^  .  -  UI 
n+ 1  n 

+  (1  -  e)(AT)([KI]^{Ul}^  +  [KIJ]^{UJ)^  +  (SI)^)  (45) 


Starting  from  guessed  UI  values,  corrections  for  UI  values  can  calculated 
using  the  Newton  iteration  algorithm: 

[J(FI)]P 

n+  1 

Iterated  solutions  can  then  be  defined  as: 


(UI) 


p+1 

n+ 1 


(UI) 


n-t- 1 


P+  1 


(47) 


Here  the  Newton  algorithm  Jacobian  is  constructed  as: 


[J(FI)] 


3{FI} 

S{UJ} 


(48) 


The  process  is  repeated  until  a  preset  convergence  criterion  is  met. 
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A  computer  program  UDCFLOW  has  been  written  using  the  above-described 
finite-element  algorithm  and  calculation  procedure  for  evaluating  distribu¬ 
tions  of  U  and  V  in  X  and  Y  domain  at  different  time  t.  To  facilitate  the 
use  of  the  program,  the  UDCFLOW  has  been  interfaced  with  the  PDA's  PATRAN 
graphics-oriented  modeling  pre-  and  post-processors. 

5.4  UDCFLOW  and  PATRAN  Interface 

UDCFLOW  is  a  finite  element  computer  program  for  simulating  transient 
flow  of  fluid  in  two  dimensional  fields.  To  facilitate  efficient  flow 
simulation,  the  UDCFLOW  has  been  intergrated  with  the  PATRAN  software  package 
for  graphic  pre-  and  post-processing  of  the  finite  element  analysis.  The 
PATRAN  program  provides  an  interactive,  graphics-oriented  environment  for 
creation  of  a  geometric  and  finite-element  model.  A  UDCPAT  interface  program, 
reads  the  PATRAN  generated  neutral  file  for  the  model  and  generates  a  data 
file  which  may  be  used  as  inputs  for  the  UDCFLOW  finite  element  computer 
program.  Simulation  results  from  UDCFLOW  runs  may  be  transferred  back  to 
PATRAN  for  post-processing  display.  As  an  example,  table  2  shows  an  example 
of  the  PATRAN  command  set  which  has  been  used  to  generate  the  cyclone  model 
shown  in  figure  14. 

TABLE  2 

AN  EXAMPLE  OF  PATRAN  INTERACTIVE  COMMANDS 
FOR  CYLONE  MODEL  GENERATION 

START: 

PATRAN 

GO 

PHASE  I— GENERATE  GEOMETRICAL  MODEL: 

GR,1,,0 
GR,4,TR,12,1 
GR, 2/3, TR, 0/10, 1/4 
LI, 1/2, ST, ,1/4, 2/3 
PA,1,2L, ,1,2 
SET, LINES, 0 
SET,CPL0T,0N 
PA, A, REV 
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(Table  2  Continued) 

PA, 1, LAB, ,1/2 

PHASE  II— GENERATE  ANALYTICAL  MOim,: 

GF, PI, G, 25/21 
CF, PA, QUAD 

DF , PI ,DIS , 3/0/0/0/0/0 , 1 

DF,P1,FORC,3/0/0/0/0/1,1,ED1T4 

SET, LABELS, OFF 

PM, 1, TAN, 1,10000,1, 0,4(0) 

PM, 2, TAN, 6(0. 025), 1000, 0.5 
PM, 3, TAN,0,0.018,2, 2,0. 005,3(0) 

PM,4, TAN, 3(5),3,2(0), 100000,30 

PHASE  III— CREATE  NEUTRAL  FIIE: 

INTERFACE 

NEUTRAL  FILE 

CREATE  OUTPUT 

ENTIRE  FILE 

OUTPUT  FILE 

OUTPUT  PHASE  I, YES 

OUTPUT  PHASE  II, YES 

NEUTRAI  FILE  PATRAN.OUT  CREATED 

STOP  PATRAN 

PHASE  IV— TRANSLATE  NEUTRAL  FILE  AND  UDCFLOW  RUN: 
RUN  PATUDC 
RUN  UDCFLOW 

PHASE  V— POST-PROOESSING  DISPLAY: 

RESTART  PATRAN 
OUTPUT  RESULTS 
USE  EXTERNAL  DATA  FILE 
VECTOR  PLOT 

USE  DATA  OTHER  THEN  INDICATED 

PLOT  RESULTANTS 

COMPONENTS  1  AND  2 

USE  UDCFLOW  OUTPUT  FILE  RESULT.DAT 

PLOT  MODEL 

PLOT  GRAPHS 
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5.  Results  and  Discussion 

The  algorithm  developed  can  theoretically  be  used  to  simulate  the  two- 
dimensional  atmospheric  vortex.  However,  because  of  the  vortex  scale  in 
nature,  supercomputer  will  be  required  to  simulate  the  atmospheric  vortices. 
The  author  has  so  far  used  only  the  UDC's  VAX-8650  computer  for  test  runs  of 
the  program  UDCFLOW.  Several  test  runs  have  been  made  for  problems  of 
academic  interests.  Runs  were  first  made  for  the  one-dimensional  fluid  flow 
and  thermal  diffusion  problems.  Several  examples  of  these  runs  have  been 
presented  and  discussed  in  a  paper  by  Chi,  Berhanu  and  Ranje  (1989)  and  a 
reprint  of  the  report  is  appended  to  the  report  (see  Appendix  B). 

A  very  popular  two-dimensional  flow  problem  is  the  so-called  driven 
cavity,  e.g.,  see  Bozeman  and  Dalton  (1973),  Kawahara  and  Okamota  (1976),  and 
Baker  (1983).  Figure  15  shows  a  sketch  of  the  driven  cavity  with  the  nondi- 
mensional  lenct'  by  height  being  one  by  one.  The  cavity  fluid  is  initially  at 
rest,  and  the  top  lid  of  the  cavity  is  started  impulsively  with  a  velocity 
which  is  used  as  the  reference  velocity  v^.  Reynolds  number  for  the  cavity 
flow  which  is  defined  as  (v  1  /v  )  is  assumed  to  be  100.  At  the  test  run,  the 
flow  domain  is  divided  into  a  mesh  of  20x20.  Figures  16  through  20  show 
snapshots  of  the  flow  field  at  the  nondimensional  time  (t)  equal  to  0.0,  0.4, 

1.0,  3.0  and  5.0,  respectively.  It  can  be  seen  in  these  figures 

that  as  time  r  advances  from  0  to  3  penetration  of  the  flow  depth  increases. 
As  time  advances  further  from  t  equal  to  3,  there  is  little  change  in  the 
cavity  flow.  That  is,  steady  state  of  the  flow  has  been  approached  at  t  equal 
to  3. 

In  conclusion,  a  general  finite  element  computer  program  UDCFLOW  has  been 
developed  to  simulate  transient  flow  develpment  in  the  two-dimensional  domain. 
The  program  has  been  integrated  with  the  PDA's  graphics-oriented  PATRAN 
software  package  for  finite-element  analyses  pre-  and  post-processing. 
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FIGURE  IS  SCHEHATIC  DIAGRAH  OF  A  DRIVEN  CAVITY 

SHOUING  INTIAL  AND  BOUNDARY  CONDITIONS 
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FIGURE  1?  DRIVEN  CAVITY  SOLUTION  AT  RE=100  AND  TAO= 
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FIGURE  Ifl  DRIVEN  CAVITY  SOLUTION  AT  RE=1DD  AND  TAO= 
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FIGURE  n  DRIVEN  CAVITY  SOLUTION  AT  RE^IUD  AND  TAO= 
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FIGURE  20  DRIVEN  CAVITY  SOLUTION  AT  RE=1D0  AND  TA0= 


Runs  of  the  program  for  several  interesting  problems  have  been  made  by 
the  author  and  his  co-workers.  Effort  is  now  being  made  by  the  author  for 
accessing  to  Super-Computers  to  test  applicability  of  the  program  to 
full-scale  tropical  cyclones  in  atmospheres. 

6.  CONCLUSIONS 

Under  the  ONR  Research  Grant  No.  N00014-84-0610 ,  UDC  has  undertaken 
a  basic  research  on  dynamics  of  marine  atmospheres  in  tropical  cyclones 
with  an  objective  of  improving  understanding  of  turbulence  in  tropical 
cyclones  and  efficiency  of  computer  simulation  of  tropical  cyclones.  For 
this  objective,  UDC  has  performed  the  following  four  phases  of  effort: 

I.  Development  of  mathematical  models  for  marine  planetary 
boundary  layers  of  tropical  cyclones  over  the  sea  surface. 

II.  Simulation  of  marine  planetary  boundary  layers  of  tropical 
cyclones  over  the  sea  surface, 

III.  Tests  in  laboratory  of  air  vortex  flows  over  the  water  surface. 

IV.  Development  of  finite  element  computer  programs  for  simulating 
marine  planetary  boundary  layers  of  tropical  cyclones. 

Through  these  phases  of  effort,  several  accomplishments  have  been 
made.  A  turbulent  theory  for  processes  of  the  transport  of  heat, 
moisture  and  momentum  in  the  marine  planetary  boundary  layer  has  been 
developed.  The  theory  uses  the  second-order  turbulent  closure  model  and 
calculates  the  turbulent  stress  components,  and  enthalpy  and  moisture 
fluxes.  In  addition,  a  simplified  k-e  turbulent  model  is  derived  from 
the  detailed  turbulent  flux  model  to  increase  the  calculation  efficiency. 
A  computer  program  based  upon  the  finite-difference  numerical  procedure 
has  been  coded  to  calculate  the  turbulent  transport  of  heat,  moisture  and 
momentum  in  the  marine  planetary  boundary  layer  of  tropical  cyclones. 

To  increase  confidence  in  turbulence  models  and  computer  simulation, 
a  vortex  chamber  with  air/water  interface  has  been  designed  and  set  up  in 
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the  UDC  laboratory.  The  chamber  has  been  instrumented  with  both  the  hot- 
film  and  laser  doppler  velocimetjer.  Tests  have  so  far  been  made  only  with 
the  hot-film  anemometer  and  with  air/solid-wall  interface,  and  data  from 
these  tests  with  both  the  smooth  and  rough  solid  walls  are  in  excellent 

agreement  with  the  computer  simulation  results. 

To  increase  simulation  efficiency,  a  major  effort  has  been  made  in 
development  of  a  graphics-oriented  interactive  finite-element  model  for 
simulating  transfer  of  heat,  moisture  and  momentum  in  vortex  flows.  A 
finite-element  analysis  computer  program  (UDCFLOW)  has  been  developed  at 
UDC.  Solutions  have  so  far  been  obtained  for  transient  one-dimensional 
heat  and  momentum  transfer  and  two-dimensional  momentun  transfer  equa¬ 
tions.  In  addition,  the  UDCFLOW  program  has  been  integrated  with  the 
PDA's  graphics-oriented  PATRAN  software  package  for  finite-element 
analysis  pre-  and  post-processing.  Successful  runs  of  the  program  for 
several  interesting  problems  have  been  made  by  the  author  and  his  co¬ 
workers.  Extension  of  runs  to  problems  of  vortex  flow  in  marine  atmos¬ 
pheres  will  need  more  extensive  computer  facility  than  the  VAX-8650  which 
the  author  has  been  using.  Effort  is  now  being  made  by  the  author  to 
access  to  "Super-Computers"  to  test  applicabiliy  of  the  program  to  study 
full-scale  tropical  cyclones  in  marine  atmospheres.  In  addition,  inves¬ 
tigations  are  being  made  of  possible  permatations  of  index  structure  for 
the  large-scale  Jacobian  '  metrices  so  as  to  decrease  the  central  memory 
reg'jirements  for  the  Jacobians  and  to  reduce  the  CPU  time  for  their 
evaluation . 
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ABSTRACT 


■  surface  total  moisture  flux 


This  paper  presents  a  numerical  simulation  of 
turbulent  vortex  bounrtary  layers  on  the  water  surface. 
The  model  of  turbulence  employed  is  one  that  the  tur¬ 
bulence  energy  and  its  dissipation  rate  are  calculated 
hy  w.iy  (if  tc.insport  i'i|ii.i  c  i  oils  wliicli  .in*  solved  siimil- 
taneouslv  with  the  vortex  boundary  layer  mo. in  flow 
equations  for  conservations  of  mass,  momentum,  energy 
and  moisture.  The  model  considered  includes  influen¬ 
ces  of  the  buoyancy,  Coriolis  and  centrifugal  forces, 
the  air/water  interface  roughness,  and  the  interface 
water  droplet  entrainment.  Although  lack  of  turbulen¬ 
ce  measurements  over  the  water  surface  prevents  a  com- 
lete  comparison  of  the  present  theory  with  experiments, 
good  agreement  is  obtained  with  available  experimental 
data  for  the  vortex  boundary  layer  over  simulated  wa¬ 
ter  surfaces. 

NOMENCLATURE 

C^  ■  constant  pressure  specific  heat 

C  'a  constant  in  equation  for  calculating  the  water 
surface  roughness 

f  •  Coriolis  parameter 


r  ■  radial  distance 
Ri  •  Richardson  number 

S  •  mi'jn  dry  s|>ucific  viithulpy  defined  as  C  IT-(T, 
-gt/Cp)l 

S  ■  virtual  dry  specific  enthalpy  defined  as  S(  !♦ 

''  t.609Q^-Q) 

T  =  temperature 

T,  “  standard  temperature 

U  •  mean  velocity  component  in  radial  di.“ction 
“  surface  frictional  velocity 

V  “  mean  velocity  component  in  tangential  direction 

V  «  mean  tangential  velocity  at  the  vortex  outside 

radius 

W  »  mean  velocity  component  in  vertical  axial  direc¬ 
tion 


g  •  gravitational  acceleration 

H  ■  mean  moist  specific  enthalpy  defined  as  (StQ^) 
H.J  “surface  moist  static  energy  flux 
k  “  turbulence  kinetic  energy 
L  •  water  latent  heat  of  vaporization 
Q  “  mean  total  moisture  mixing  ratio  defined  as 

Qj  “  liquid  water  mixing  ratio 
0^  “  water  vapor  mixing  ratio 


X  “  dependent  variables 

z  ■  vertical  axial  distance 

z,  “  water  surface  roughness  length 

X,  ■  first  layer  of  air  next  to  the  air/water  inter¬ 

face 

a  “a  coefficient  in  finite  difference  equations 

8  “  expansion  coefficient  or  a  coefficient  in  finite 

difference  equations 

9  “a  function  of  Richardson  number  defined  in  text 
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Tf  “  •  quantity  defined  aa  (L/C  )(3C}  /3T  )  or  a  coef¬ 
ficient  in  finite  differen§e  eqSatiSna 

4  “a  coefficient  in  finite  difference  equations 

A  •  boundary  layer  thickness 

c  •  turbulence  energy  dissipation  rate 

t  ~  a  quantity  defined  as  C  T,/L 

P 

4  ■  azimuth  angle 

•  turbulent  kinematic  viscosity 
0  •  Prand t 1 /Schmid t  number  for  turbulent  diffusion 

T  'a  function  of  Richardson  number  defined  in  text 
SUBSCRIPTS 

g  >  saturated  vapor 

h  “  pertaining  to  diffuaional  transport  of  moist 
specific  enthalpy 

k  •  pertaining  to  diffusional  transport  of  turbu¬ 
lence  energy 

m  •  values  of  the  radial  position  where  the  main 
vortex  tangential  velocity  is  maximum  or  the 
vertical  grid  position  in  a  finite  difference 
scheme 

H  "  grid  position  at  the  vortex  top  boundary 

n  •  the  horizontal  grid  position  in  a  finite  dif¬ 

ference  scheme 

N  “  grid  position  at  the  vortex  outside  radius 

q  •  pertaining  to  diffuaional  transport  of  moisture 

s  ■■  pertaining  to  diffusional  transport  of  dry  spe¬ 

cific  enthalpy 

X  •  pertaining  to  the  dependent  variable  X 

c  "  pertaining  to  diffuaional  transport  of  turbu¬ 

lence  dissipation  rate 

*  ~  lower  boundary  of  an  air  vortex  at  the  nir/water 

interface 

'  ”  air  layer  next  to  the  air/water  interface 

•  "  outside  radial  boundary  of  an  air  vortex 

♦  ■  top  boundary  of  an  air  vortex 
INTRODUCTION 

In  their  earlier  papers,  the  present  author  and 
his  co-workers  reported  several  numerical  simulation 
results  for  the  turbulent  vortex  boundary  layer  on 
smooth  solid  surfaces,  using  the  Karman's  integral 
method  Ml.  The  mixing  length  theory  12),  and  the 


'  The  number  in  brackets  refer  to  the  references  at 
the  end  of  the  paper 


turbulent  energy  method  (31,  respectively.  Recent  ad¬ 
vances  in  turbulence  modeling  and  numerical  methods 
(e.g.,  sec  Harlow  and  HirtlA],  Gibson  and  Spalding  (51, 
Hanjalic  and  Launder  (61,  Mellov  and  Vamada  (71,  Le- 
wellen  (8),  and  Chi  l9))  have  made  it  a  common  prac¬ 
tice  to  solve  simultaneously  transport  equations  for 
turbulence  quantities  with  conservation  equations  for 
mean  flows.  Many  parameters  for  complex  turbulent 
flows  can  now  be  simulated  by  numerical  methods.  In 
this  paper,  the  author  presents  a  numerical  simulation 
of  turbulent  vortex  boundary  layer  flows  over  the  wa¬ 
ter  surface  like  that  shown  in  Fig.  I. 


Fig.  1  Sketch  of  the  Vortex  Flow  Over  a  Water  Surface 

The  flow  situation  shown  in  Fig.  1  occurs  in  na¬ 
tural  ati.ospheres  such  as  the  marine  planetary  boun¬ 
dary  layer  of  tropical  cyclones  as  well  as  engineer¬ 
ing  devices  such  as  wet  cyclones  and  thermal  proces¬ 
sors.  It  presents  many  challenging  problems  which  are 
of  interests  to  researchers  in  areas  of  heat  and  mass 
transfer,  fluid  dynamics  and  marine  meteorology.  The¬ 
se  problems  include  the  boundary  l.ayer  lic.it,  moisture 
and  momentum  transfer;  turbulence  production,  dissi¬ 
pation  and  diffusion;  strong  influence  of  the  gravita¬ 
tional,  Coriolis  and  centrifugal  forces;  and  influence 
of  air-water  interface  roughness  and  liquid  droplet 
entr.aiiunent . 

The  model  of  turbulence  employed  in  this  paper  is 
one  which  calculates  the  turbulent  kinetic  energy  k 
and  its  dissipation  rate  c  by  way  of  their  transport 
equations.  The  turbulent  kinematic  viscosity  will  be 
calculated  from  k  and  c  values  (lOj, 

v^  •  0.09k*/c  (1) 

Mean  conservation  equations  which  will  be  solved  si¬ 
multaneously  are  those  not  only  for  balances  of  mass  and 
momentum  but  also  for  balances  of  the  moist  static 
energy  and  the  total  moisture  mixing  r.itio  so  as  to 
en.iblc  the  model  to  simul.atc  both  the  air/w.iter  inter¬ 
face  evaporation  and  the  interface  water  droplet  en¬ 
trainment.  The  Coriolis  force  due  to  earth  rotation 
and  the  centrifugal  force  due  to  the  flow  circulation 
are  included  in  formulating  the  model  momentum  balance. 
Influences  of  buoy.ancy  force  .ind  surf.ice  roughness  are 
accounted  for  in  formilating  the  model  turbulence 
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CransporCac ion  and  the  model  flow  boundary  conditions. 

In  subsequent  two  sections  of  this  rei>ort  ,  gover¬ 
ning  equations  for  the  model  are  derived  along  with 
their  relevant  initial  and  boundary  conditions  and 
th«'ir  s<^)luci«.>n  procedure,  .Mtd  results  ure  presented 
and  d i scussed . 


THE  MODEL 

Covcrnifi>;  Etiuations 

Coord i nates  are  chosen  so  that  z  axis  is  vertical, 
r  axis  is  in  the  r.idinl  direction  and  i  is  tlic  aeiinuiii 
angle.  When  tlie  Boussinesq  and  boundary  layer  appro¬ 
ximations  and  the  axisymmetrical  and  quas i -s teady  f low 
assumptions  are  used,  governing  eqtiationa  for  conser¬ 
vations  of  mean  mass,  mean  motion  rad i a  1 -d i rec t ion  mo¬ 
mentum,  mean  motion  tangent ia I -d i rec t ion  momentum, 
mean  motion  moist  static  energy,  mean  motion  total  mo¬ 
isture  mixing  ratio,  mean  turbulence  kinetic  energy 
and  mean  turbulence  dissipation  rate  become 
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(3) 

iU) 

(5) 

(6) 

(7) 
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In  the  above  equation! ,  U,  V  and  W  are  the  mean 
velocity  components  in  the  r,  e  and  z  directions;  and 
H  is  the  mean  moist  specific  enthalpy,  Q  the  mean  total 
moiature  mixing  ratio,  k  the  mean  turbulence  kinetic 
energy  and  c  the  mean  turbulence  energy  dissipation 
rate.  The  turbulent  viscoaity  appearing  in  die  diffu¬ 
sion  term  is  calculated  by  equation  (1).  Tlie  turbulent 
exchange  coefficients  for  diffusions  of  the  moist  sta¬ 
tic  energy,  the  total  moisture  ratio,  the  virtual  dry 
spefiric  enthalpy,  the  turbulence  kinetic  energy  and 
the  turbulence  dissipaciun  rate  are  related  tu  die  tur¬ 
bulent  viscosity  through  the  turbulent  Prandt l-Schmidt 
numbers  wliich  are  all  tuppoaed  to  take  on  the  follow¬ 
ing  constant  values: 


predicted  X  distributions.  However,  the  buoyancy  terms 
in  equations  (7)  and  (8)  contain  a  derivative  of  S 
which  is  not  one  of  the  directly  predicted  variables. 

In  order  to  relate  (oS  /3z)  to  the  predicted  variables 
H  and  Q,  we  follow  a  derivation  discussed  by  Hoeng  and 
Ar.'ikawa  [II],  which  yields  the  following  re  1  a  t  ionsh  i  ps  : 


3S 

V 

17 

as 


,  3H 

II 


♦  -  l)L3Qj  region)(IO) 


a  3  z 
‘I 


V  ,  (  I  ♦  1 . 609tT )sH  cLSQi  ,,  ,  .  .  , 

- — “0  ! — 7-, - ; - (for  cloudy  region)  (11) 

3z  8  (1*1)0,32  o3z 

h  <1 


Here  y  is  equal  to  (L/C  )(3Q  /3T  ),  t  is  equal  to  C  T  / 
P  g  g  P  o 

L,  and  sufix  g  stands  fur  the  saturated  water ennd i t i on. 
Equations  (2)  through  (8)  together  with  the  auxiliary 
functions  discussed  above  now  form  a  model  for  turbu¬ 
lent  vortex  flows  over  the  water  surface.  In  addition, 
the  boundary  layer  model  requires  specification  of  ini¬ 
tial  conditions  and  boundary  values  at  the  edges  of  the 
flow. 


Bottom  Boundary  Conditions 

At  the  air/water  interface,  both  the  water  surface 
roughness  and  Che  surface  layrt  stability  influence  the 
air  vortex  bottom  boundary  conditions.  The  water  sur¬ 
face  roughness  length  can  be  calculated  by  the  equation 
112): 

2  •  C  Ui/g  (12) 

where  U.^,  is  the  surface  frictional  velocity,  g  Che  gra¬ 
vitational  acceleration,  C  s  constant  whose  value  de- 
pends  on  the  extent  of  the  flow  field.  Tor  example, 
for  the  air/sea  interface  it  is  equal  to  1.6x10"’,  for 
Che  laboratory  wave  channel  [12]  it  is  equal  to  1.l2x 
10"’,  and  for  the  UDC  vortex  chamber  [9]  which  will  be 
discussed  later  in  Che  paper  it  is  equal  to  5.58x10"’. 

For  numerical  modeling,  the  boundary  layer  in  the 
vortex  flow  is  divided  into  horizontal  layers  and  con¬ 
centric  cylindrical  surfaces.  The  first  bottom  layer 
will  be  at  a  height  z*z,.  U,  V,  H  and  Q  values  at  the 
surface  layer  can  be  calculated  by  tlie  equations  [  13,  IA  1: 


(U|*  V’)^.  2 

.50U  [ln(^-A)  ♦  T(Ri,  )1 

*  • 

(13) 

(H,-  H.)  •  1 

.85H.  [  In(^-L)  ♦  e(Ri ,  )  1 

(  IA) 

(Qi-  Q.)  ■  1 

.85Q,^1  In(^-k)  ♦  e(Ri , )  1 

2  • 

(  15) 

Here  U^,  and  are  the  surface  frictional  velocity, 

the  surface  moist  static  cnerj;/  flux,  and  the  surface 
total  moisture  flux,  respectively;  uuJ  the  Uiclurdsou 
number  Ri  is  defined  js  follows: 

V 

Ri  -  i0.4y(- 
0 

i)(^),l/(uis^; 

(  16) 

■ 


(o,.,  0  ,  o  ,  0,.,  o  )  -  (0.9, 0.9, 0.9, 1.0, 1.3)  (9) 

h  q  a  k  c 

In  Che  above  equation  set,  the  turbulence  stresses 
snd  fluxes  have  been  modeled  by  the  expression  [(v  /o  ) 
(3X/32)|  in  which  X  represents  any  mean  quantityof^a  * 
dependent  variable  value.  Values  of  [(v  /o  )(3X/32)I 
in  equations  (2)  through  (8)  can  be  calculated  from  Che 


The  stability  functions  T(Ri)  and  6(Ri)  are  calculated 
by  the  eqiiat  ions  : 

T(Ri )  -  A  .  7Ri  [ for  Ri  >  Ij 

-  -ln(0.5[l*(l-l5Ri)'*l)-21n(0.5[l*(l-t5Ri)^)) 

♦  2tan"'(  l-l5Ri)^-0.5ll  I  for  Ri  <  1)  (W) 


629 


e(Ri)  -  6.35Ri 


I  for  Ki  > 


-  -21nt0.5ll*(l-9IU)  Ml  ICor  Ri  ^  11  (I 

k  and  c  values  at  the  botton  boundary  layer  can  be 
calculated  by  the  equations  1151: 

k  -  3.33U'(l*i.7Ri,)*’  Ifor  Ri  ^  1] 

-  3.33Ui( l-l5Ri,)^  Ifor  Ri  ^  1 1  (1 

c  -  0.4I2k^^^z7'  (3 


Top  Boundarv  Conditions 


H  -  H,  -  1.8511.  lln(  -)  ♦  0(Ki)) 


Q  -  Q,  •  l.85Q.|ln(— )  *  0(Ri)) 


Here  Ri,  T  and  0  values  are  deterniined  by  ef;uations 
(16)  through  M8),  respectively.  The  side  turbulence 
energy  profile  may  be  represented  bv  the  polynomial 

n9l: 

k  -  ((k,)  /[I  -  2(^L)i.  3(lL)>])ll-2(v)'*3(i)’l  (31) 

•  A  A  A  A 


Mean  velocity  components  for  the  m.iin  vortex 
flow  at  and  all  r  hnve  been  calcul^ited  by  numo'^ 

roua  workers,  e.g.,  see  RootllCl,  Sullivan  ll?!  and 
Lewelleo  [18).  The  result  can  be  represented  accura¬ 
tely  by  tlie  equations: 


(U  )  (r  ) 


1 .  ^  V  r 

-  - -JUB  (I  -  exFl-l.253(^)’)) 


in  which  (k,  )  is  tlie  bottom  boundary  turbulence  energy 
value  at  r*r  It  can  be  calculated  by  e<juation  (19). 
V.ilue  for  the  side  bouf»dnry  turbulence  dissipation  rate 
is  related  to  the  turbulence  energy  value  by  the  equa- 
t  Ion : 

.  .  .  ,  , ,,.3/2-1  (32) 


Solution  Procedure 

The  conservation  equations  (3)  through  (S)  for  U, 
V,  H,  Q,  k  and  c  have  a  parabolic  form  of 


where  suffix  *  represents  values  at  the  top  boundary. 
Suffix  •  represents  values  at  large  radius,  and  suffix 
m  represents  values  at  a  position  in  the  main  flow 
where  tlie  vortex  tangential  velocity  is  maximum.  In 
the  main  flow,  the  mean  moist  static  energy  H  and  the 
mean  total  moisture  ratio  Q  values  are  expected  to  be 
linear  with  respect  t<>  7. .  Hence  tlieir  top  bound.iry 
conditi  .8  are  as  follows: 


=  .  Source  (33) 

3r  3z  3z  0  Sz 
X 

Severjl  finite-difference  methods  of  solution  for  equa¬ 
tion  (33)  are  available  in  the  literature  [201.  We  ha¬ 
ve  written  n  comi'iifer  protir.nm  followim;  on  extension  >'( 
tin-  implicit  procedure  civiui  hy  c’r.ink  .ind  N  i ii,' 1  stni  I.' I  I  . 
Using  a  rectangular  grid  system  on  the  (r,z)  plane  with 


The  values  of  k  and  c  at  the  top  boundary  can  be  de¬ 
termined  from  the  degenerative  forms  of  equations  (7) 
and  (8)  which  result  when  gradients  with  respect  tc  z 
are  set  to  zero.  That  is, 


r  ♦  4r 
n  n 


z  ♦  4z 
m  n 


(n»l  ,2 . N) 

(m*  1.2 . >1) 


The  equation  set  (33)  for  different  quantities  of  the 
model  independent  variables  are  approximated  at  the  po¬ 
int  (r  .,,z  ) and  written  in  the  form' 
n*  1  m 


n  n  .n„n 

o  X  B  X 

m  m— I  m  m 


m  ra-  I 


(U  )(k  ■) 


Side  Initial  Values 

For  the  side  boundary  initial  values,  the  boun¬ 
dary  thickness  B,  the  radial  velocity  amplitude  value 
E  and  the  surface  friction.rl  velocity  U  at  the  outside 
radius  r  •  r„  may  first  be  calculated  by  a  Karmaii's 
integral  method  |l|.  Initial  values  for  U,  V,  II  and 
Q  St  r*r.  and  all  z  can  then  be  repreaenCed  by  the 
(ol lowing  functiona  (  I ,  lA  |  : 


V  - - r(ln(-=-)  T(Ri)I 

(HE’)' 

(U  ) 

U  -  EV(S[1  ♦  coa(:^)|)(l  -  211  -  Etv-)"  l<f) 

(U .) 


*  II  - 


Here  at  n*N,  values  of  X  are  equal  to  the  side  initial 
values  at  r"r_;  .-t  m"  I ,  values  of  X  are  equal  to  the 
bottom  boundary  values  at  z”z,;  and  at  rn'M,  values  of  M 
are  equal  to  the  top  boundary  values  at  Integra- 

l.inr.  lupi.'i  (  i  vui  (II)  I  ri,m  n-N-  I  with  m  I  litm  I  lu  M,  v.i- 
mes  of  x'”  j  with  m  from  I  to  M  c.in  be  di'tormi  ned  .  V.i- 
luc  of  II  IB  then  reduced  by  uiu*  tu  dutiTini  iiu  .\  Iho 

process  is  continued,  until  it  covers  the  whole‘'(r,z) 
dop.:in  for  the  problem.  Structural  details  of  the  pro¬ 
gram  can  be  found  in  a  paper  by  Chi  and  Glowacki  [2|. 

RESULTS  AND  CONCLUSION 

-  In  order  to  ettabliah  coi.fidence  in  Che  present  mo¬ 
del,  tome  measurementa  of  the  turbulent  vortex  flow  ha¬ 
ve  recently  been  made  (on  a  UDC  vortex  teat  chamber)  by 
the  author  and  a  co-worker.  Details  of  Che  experimental 
set-up  (Fig.  2)  were  given  by  Chi  and  Hinds  [91.  Brief¬ 
ly.  the  test  chamber's  overall  dimension  is  55-cm  diame¬ 
ter  by  90-cm  high.  It  generates  an  air  vortex  of  30.^8- 
cm  diameter  by  15-cm  high.  The  vortex  is  generated  bv 
forcing  air  through  2A  evenly  sp.iccd  vanes  of  15-cm  high 
placed  on  a  30.A8-cm  diameter  pitch  circle.  The  swirl¬ 
ing  air  disch.arges  from  the  ch.imbcr  through  a  3.81-cm 
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diamecer  liole  located  at  center  of  the  chamber'a  top 
disk.  Bottom  of  tlie  swirling  air  is  in  contact  with 
s  pool  of  water  which  can  be  maintained  at  a  constant 
depth  and  at  desired  temperatures.  A  scaled  needle 
probe  may  be  inserted  vertically  from  the  bottom  into 
the  chamber  at  different  radial  positions  to  measure 
mechanically  the  wave  heights. 


AIR  OJTLET 


AIR  INLET 


AIR  FLOW 


AIK  FLOW 
GUIDE  VANES 

’‘JTFUFACE 

AiR  FLOW 

TURBULENCE 

PROut, 


Fig.  2  Schematic  of  a  Laboratory  Vortex 
Chamber  with  Air/Water  Interface 

Because  the  chamber  is  not  currently  instrumented 
to  compensate  for  influences  of  water  entrainment  on 
the  air  flow  raeasurenenc,  direct  measurements  of  air 
velocity  over  the  actual  water  surface  have  not  been 
maoe.  Solid  surfaces  of  different  roughness  heights 
were  used  to  simulate  physically  the  air/water  inter¬ 
face  roughness.  The  hori7.onc.1l  component  of  velocity 
and  its  direction  in  the  vicinity  of  Che  rough  solid 
anrfjce  have  been  measured  by  a  directionally  aenai- 
tive  wedge-ahaped  hot-film  probe.  The  hot  film  waa  on 
the  axis  of  Che  cylindrical  stem  of  Che  probe,  and  Che 
probe  was  tr.Tvrrsed  .ixi.nlly  thrinit;h  the  bottom  plate 
.It  different  r<id  i  i  so  th.it  the  h«>ri40ut.il  vclocit.i%*s 
at  different  radii  and  heights  could  be  measured.  Fi¬ 
gures  3  and  U  show  tb  measured  vortex  boundary  layer 
velocity  profiles,  U  and  V,  respectively,  versus  z  at 
several  radii  and  with  the  surface  roughness  values 
indicated . 


RADIAL  INWARD  VELOCITY,  -U 

Fig  3  Predicted  and  Measured  Radial  Velocity  on  Solid 

Surface.  (Keys:  o  &  - ,  Measured  S  Calculated 

values  on  Smooch  Surface;  o  &  Measured  U 

Calculated  Values  on  Rough  Surface  of  Sand  Crain 
Mesh  fll20) 


0  0  0.0  0.0  0.0 

TANGENTIAL  VELOCITY,  V 


Fig.  4  Predicted  and  Measured  Tangential  Velocity  on 

Solid  Surface.  (Keys:  o  &  -  ,  Measured  &  Cal¬ 
culated  Values  on  Smooth  Surface;  o  &  Mea¬ 

sured  &  Calculated  Values  on  Rough  Surface  of 
Sand  Crain  Mesh  i?120) 

As  numerical  examples,  the  computer  program  for 
the  present  model  has  been  run  to  simulate  the  turbu¬ 
lent  vortex  boundary  layer  flows  under  the  same  main 
vortex  conditions  as  in  Che  experiments  out  with  seve¬ 
ral  different  bottom  surface  conditions  as  follows: 

1.  Smooth  solid  surface, 

2.  Sund-roughnoss  solid  surface,  and 

3.  Actual  water  surface. 

For  amocth  solid  surf. ice  c.i  Icii  1  jC  ions  ,  values 
were  calculated  from  4,  eipul  to  U.lllv/U*,  and  for  sand 
roughness-surface  calcula  ions,  values  were  calcula¬ 
ted  from  Z|  equal  to  l/3(ith  of  the  sand  grain  height. 
Both  of  these  relationships  have  been  obtained  from 
Schlichting  l22i.  The  \>r<‘dicted  rodiul  velocity  U  and 
t.ingential  velocity  V  over  the  smooth  and  rough  surfa* 
ces  are  plotted  versus  z  for  several  radii  in  Figs.  3  and 
4  for  comparison  with  the  experiments.  The  ability  nf 
the  present  theory  to  correlate  the  measured  mean  velo¬ 
city  profiles  in  turbulent  vortex  boundary  layers  over 
smooch  and  rough  solid  surfaces  is  demonstrated  in  Figs. 
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TACENTIAL  VELOCITY 


Fig.  6  Predicted  Tangential  Velocity  Over  a  Water  Sur¬ 
face 

Mnny  intereating  clinrac  t  er  i  a  t  i  ca  for  tlic  vortex 
boundary  layer  flow  can  be  observed  from  llie  above  pre¬ 
sented  numerical  aimulations.  Near  the  a i r / wa t e r  in t er- 
face,  retardation  of  the  tangential  velocity  is  seen  in 
'’iga.  3  and  6.  Tliis  retardation  ia  accompanied  by  a  re¬ 
duction  in  the  centrifugal  force;  the  balance  between 
iresaiire  and  the  centrifugal  force  is  t  he  reby  des  t  royed  . 
Hie  flow  in  tliia  region  ia  thus  characterized  bv  the 
entrainment  of  air  from  the  main  vortex  into  tlie  boun- 


0.0  0.0  0.0  o.n 

TURBULENCE  KINETIC  ENERGY,  k 

Fig.  9  Predicted  Turbulent  Kinetic  Energy  Over  .a  Water 
Sur  f ac  e 

dary  layer  as  indicated  by  an  excess  of  Che  radial  velo 
city  inside  the  boundary  layer  over  its  value  in  the 
main  vortex  at  the  same  radius.  Effects  of  the  solid 
surface  roughness  and  the  water  surface  waviness  can  be 
seen  in  Figs.  3,  A,  5  and  6  to  increase  the  boundarv 
layer  thickness  but  to  decrease  the  abruptness  in  dece¬ 
leration  of  the  tangential  flow  and  in  acceleration  of 


Che  radial  flow.  In  addition,  the  large  radial  inward 
veloci:;'  inside  Che  vortex  boundary  layer  carries  cir¬ 
culation,  static  enthalpy,  noiscure  content  and  turbu¬ 
lence  energy  frore  large  radii  to  a  position  near  the 
axis  o:  avTtmetry  as  can  be  seen  at  snail  radius  in  Fi¬ 
gures  -  and  0  of  the  large  V,  in  Fig.  7  of  the  large 
in  Fig.  3  of  large  and  in  Fig.  9  of 

tlie  large  k. 

In  conclusion,  a  numerical  simulation  of  turbu¬ 
lent  boundary  layer  flows  on  Che  water  surface  has 
been  developed.  The  model  of  turbulence  employed  is 
one  that  the  turbulence  energy  and  its  dissipation 
rate  are  calculated  by  way  of  transport  equations 
which  ate  solved  simultaneously  with  the  vortex  boun¬ 
dary  layer  conservations  of  mass,  momentum,  moist 
static  enthalpy  and  total  moisture  ratio  for  Che  mean 
motion.  As  far  as  possible,  predicted  flows  are  com¬ 
pared  with  available  measured  mean  values  in  the  vor¬ 
tex  boundary  layer.  The  ability  of  Che  present  model 
to  correlate  the  empirical  data  is  shown  to  be  good 
(see  Figs.  3andH).  In  addition  several  interesting 
characteristics  of  turbulent  vortex  boundary  layer 
flows  over  the  water  surface  are  predicted  by  the  mo¬ 
del.  Predictions  of  Che  radial  inward  convergence  of 
moss,  circulation,  enthalpy,  moisture  and  turbulence 
energy  from  large  radii  to  a  position  near  the  axis  of 
symmetry  are  in  agreement  with  the  known  physical 
laws  . 

ACKNOT.iTIihENT 

This  work  was  supported  by  Che  Office  of  Naval 
Research.  The  author  wishes  to  express  his  appreciat¬ 
ion  to  the  Office. 

REFIRENCii 

1.  Chi,  S.  ,  Ying,  S.  I.  and  Chang,  C.  C.  , 

"The  Crpur.d  Turbulent  Soundary  Layer  of  a  Stationary 
Tornado-L;  ke  Vortex,"  TtLL'JS,  Vo  1 .  2  I  ,  1969. 

2.  Chi,  S.  W.  ,  and  Cldwacki,  C'.J.  ."Applicabili¬ 
ty  of  Mixing  Length  Tneory  to  Turbulent  Boundary  La¬ 
yers  Senesih  Intense  Vortices,"  J.  Appl.  Mech.,  Vol. 

LI,  191-. 

3.  Cni,  J.,  "'I'umerica  1  Analysis  of  Turbulent  End- 
Wall  Boundary  Layers  o:  Intense  Vortices,"  J.  Fluid 
Keen.,  Vo  1 .  8  2  ,  19  7  7  . 

L.  Harlow,  F.  H,  and  Hirt,  C.  W. ,  "Generalized 
Transport  Equations  of  Anisotropic  Turbulence,"  Los 
Alamos  Sci.  Lab.,  Rent.  No.  LA-4086,  1969. 

5.  Cibson,  M.  M.  and  Spalding,  0.  B.  ,  "A  Two- 
Equation.  Model  of  Turbulence  Applied  to  the  Prediction 
of  Heat  and  Mass  Transfer  in  Wall  Boundary  Layers," 

ASMS  Paper  No.  72-HT-I5,  1972. 

6.  Hanjalic,  K.  and  Launder,  B.  E.,  "Fully  Deve¬ 
loped  Asynmetric  Flow  in  a  Plane  Channel,"  J.  Fluid 
Meth.  ,  Voi .  5  2  ,  197  2  . 

7.  Mellor,  L.  C.  and  Yamada ,  T. ,  "A  Hierachy  of 
Turbulent  Closure  Models  for  IMjiieCary  Boundary  La¬ 
yer  s  , "  J .  Atmos .  Sci.,  Vol ,  3  1  ,  1974. 

8.  Lewellen,  W.  S.,  "Modeling  the  I.owesC  I-kn  of 
the  Acmos‘cere,"  North  Acl.^ntic  Tre.ity  Ori‘,.tn  i  za  C  i  on  , 
K.-pt  .  N...  ,'...-.’07,  lu.'ll. 

9.  u'.  i  ,  J.  and  Hinds,  C.  H.,  "A  Turbulence  Theory 
for  Vortex  Marine  Planetary  Boundary  Layers  and  Compa¬ 
rison  wit,'  Experiments,"  University  of  the  District  of 
Co'.umbia,  .-ept.  No.  HHEG-9511,  1985. 


’The  auftr's  name  has  since  1  97  7  been  changed  to  J. 
Ch  i 


10.  Jones,  'W,  P.  and  Launder,  B.  E.,  "The  Calcu¬ 
lation  of  Low-Seynolds-Number  Phenomena  with  a  Two- 
EquJCion  .Model  of  Turbulence,"  ASME  Paper  No.  72-HT-20. 
1972. 

11.  Moeng,  C.  H.  and  Arakawa,  A.,  "A  Numerical 
Study  of  a  Marine  Subtropical  Stratus  Layer  and  Its 
Stability,"  J.  Atmos.  Sci.,  Vol.  37,  1980. 

12.  Wu ,  J.,  "Laboratory  Studies  of  Wind-Wave  in¬ 
teractions,"  J.  Fluid  Mech.,  Vol.  34,  1968. 

13.  Paulson,  C.  A.,  "The  Mathematical  Repreaenta- 
tion  of  Wind  Speed  and  Temperature  Profiles  in  the  Un¬ 
stable  Atmospheric  Surface  Layer,"  J.  Appl.  Meteor., 

1970. 

14.  Deardorff,  J.  W.  ,  "Parameterization  of  the 
Planetary  Boundary  Layer  for  Use  in  General  Circula¬ 
tion  .Models,"  Mon.  Wea.  Rev.,  Vol.  100,  1972. 

15.  Businger,  J.  A.,  Wynaard,  J.  C. ,  Izumi,  Y. 
and  Bradley,  E.  F.  ,  "Flux-Profile  Relationship!  in  the 
Atmospheric  Surface  Layer,"  J.  Atmos.  Sci.,  Vol.  28, 

1971. 

16.  Rott,  N. ,  "On  the  Viscous  Core  of  a  Line  Vor¬ 
tex."  ZA.MM,  Vol.  10,  1979. 

17.  Sullivan,  R.  D. ,  "A  Two-Cell  Vortex  Solution 
of  the  Navier-Scokes  Equations,"  AIAA  Journal,  Vol.  26, 
1959. 

18.  Lewellen,  W.  S.,  "A  Solution  for  Three-Dimen¬ 
sional  Vortex  Flow  with  Strang  Circulation,"  J.  Fluid 
Mech. ,  Vol.  14,  1962. 

19.  Chi,  S.  W.  and  Chang,  C.  C. ,  "Effective  Vis¬ 
cosity  in  a  Turbulent  Boundaly  Layer,"  AIAA  Journal, 
Vol.  7',  1969. 

20.  Ames,  W.  F.  ,  "Nonlinear  Partial  Differential 
Equations  in  Engineering,"  Academic  Press,  New  York, 
NY,  1965. 

21.  Crank,  J.  and  Nicholson,  P.,  "A  Practical  Me¬ 
thod  for  Numerical  Evaluation  of  Solutions  of  Partial 
Differential  Equations  for  the  Heat  Conductions  Type," 
Proc.  Caoo.  Phil.  Soc.,  Vol.  43,  1947. 

22.  Schlichting,  H.  ,  "Boundary  Layer  Theory,"  Mc¬ 
Graw-Hill,  New  York,  NY,  1960. 


633 


APPENDIX  B- -REPRINT  OF  A  1989  PAPER 


Chi,  J.,  E.  Berhanu  and  M.  H.  Ranje,  1989:  A  Finite 
Element  Computer  Program  for  Dynamic  Simulate  of 
Thermal  and  Fluid  Systems,  Proc.  1989  ASEE  Annual 
Conference,  2,  788-791. 


B-1 


Session  2633 


A  Finite  Element  Computer  Program  for  Dynamic 
Simulation  of  Thermal  and  Fluid  Systems 


Joseph  Chi 

Eskinder  Berhanu 

Mohammad  H.  Ranje 

Mechanical  Engineering 

University  of  The  District  of  Columbia 

Washington,  DC 


ABSTRACT 


Computer-aided  modeling  of  thermal  and  fluid 
systems  has  gained  acceptance  both  in  scientific 
research  and  industrial  design.  Typical  CAD/CAE 
systems  for  thermal/fluid  modeling  include  a 
finite  element  analysis  code  and  graphic  facili¬ 
ties  fcr  pre-  and  pcst-processing .  The  authors 
have  developed  at  UDC  several  thermial/f lu ia  sys- 
terr.s  simulation  computer  programs  [1-3).  A  fi¬ 
nite  eie.me.-.t  analysis  co.mputer  program  UDCFLOW 
has  now  been  integrated  with  the  PATRA-N  graphic 
software  package  for  finite  element  pre-  and 
post-processing.  In  order  to  introduce  students 
to  the  new  material  and  to  demonstrate  to  stu¬ 
dents  that  the  engineer  must  stay  abreast  of 
current  developments,  'JL'CFLCW  and  PATfuAk  software 
packages  are  being  used  to  teach  students  at  UDC 
the  transient  characteristics  of  fluid  flow  and 
heat  transfer. 


INTRODUCTION 

Traditional  introductory  fluid  mechanics  and 
heat  transfer  courses  cover  simple  fluid  flow  and 
heat  transfer  problems  in  a  seg-uential  and  com- 
partmental  fashion.  Studies  of  the  transient 
fluid  flow  and  heat  transfer  problems  are  often 
limited  to  those  problems  for  which  either  the 
closed-form  solutions  or  tabulated  functions  are 
readily  available.  By  providing  students  with  a 
finite  element  f luid/therrr,al  analysis  computer 
program  integrated  with  the  graphic  pre-  and 
post-processor,  the  depth  and  breadth,  as  well  as 
the  students'  mastery  of  m.aterials  covered  in  the 
fluid/thermal  courses,  can  be  increased  over  the 
traditional  approach. 


permitting  them  the  routine  solution  of  more 
involved  problems,  and  developing  their  ability 
to  modify  programs  for  dealing  with  new  phe- 
or.-na.  For  this  objective,  the  authors  have 
developed  a  UDCFLOW  finite  element  fluid/thermal 
analysis  computer  program. 

In  this  paper,  the  UDCFLOW  computer  progra-t  is 
described  first.  Interface  of  UDCFLOW  with 
PATRAN  finite  element  graphic  pre-  and  post¬ 
processing  software  package  is  then  reviewed. 
Finally,  examples  of  coursework  that  uses  the 
UDCFUDW/PATRAN  softwares  are  presented. 

UDCFLOW  FINITE  ELEMENT  COMPUTER  PROGRAM 
AND  PATUDC  FOR  FATRAN/UDCFLOW  INTERFACING 

UDCFLOW  is  a  general-purpose  computer  progra-m 
which  is  being  developed  at  the  University  of  the 
District  of  Columbia.  To  facilitate  efficient 
simulation  under  a  wide  variety  of  conditions, 
the  UDCFLOW  has  been  intergrated  with  the  PATRAN 
software  package  for  graphic  pre-  and  post-pro¬ 
cessing  of  the  finite  element  analysis.  Although 
the  UDCFLOW  computer  program  is  being  developed 
for  simulating  3-D  heat,  mass  and  momentum  trans¬ 
fer.  For  use  in  the  undergraduate  classroom,  a 
version  of  UDCFLOW  for  solving  the  general  1-D 
heat,  mass  and  momentum  transfer  problems  has 
been  employed. 

The  governing  equation  for  the  general  1-D 
convection  and  diffusion  problem  can  be  written 
as: 


^lI.ftuiJ..i^(T-T)  (I) 

at  ax  pC  e  pc  ax  pC  ^x 
p  p  p 


Generalised  computer  programs  which  allow  a 
wide  variety  of  fluid/thermal  systems  to  be  mo¬ 
deled  and  simulated  have  recently  been  developed 
and  they  are  becoming  commercially  available. 
Examples  of  these  programs  include  Fidap, 
Flotran,  Fluent  and  Phoenics  (4,5).  However, 
source  codes  for  these  program, s  are  not  available 
to  students,  and  no  single  software  is  always 
superior  to  the  others.  For  classroom  adaptation, 
it  is  necessary  to  have  a  program  developed  by 
the  instructor  for  the  instructional  purpose  of 
relieving  the  students  of  arithmatical  drudgery. 

Professor  and  Cnairrr,an 
'  Teaching  .Assistants 


and  the  required  boundary  conditions  for  the 
problem  are  as  follows: 

At  x=0:  -I(kA/pC  )5T/5x)  =  y  T  ♦  o  (2) 

p  O  O  6  f  O  O 

At  x=L:  l(kA/pCp)3T/8xl^=  (3) 

In  the  above  equations,  t  is  the  time,  x  the 
coordinate.  T  the  temperature,  U  the  velocity,  A 
the  cross-sectional  area,  A  the  perimeter,  c  the 
density.  Cp  the  specific  heat,  k  the  thermal 
conductivity,  h  the  heat  transfer  coefficient,  q 
the  heat  source  per  unit  volume  and  T  the  envi- 
ronmental  temperature,  and  y^,,  c^,  and  are 
constant  values  at  the  domain  boundary. 
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The  (jalerkir.  finite  element  spatial  discreti¬ 
zation  of  equation?  (1)  to  (3)  lead?  to  a  system 
of  coupled  nonlinear  ordinary  differential  equa¬ 
tions  : 

ICHc}'  n  |Kl{u)  •  {S}  =  0  (4) 

Here,  [C]  is  the  capacitance  matrix,  |K)  the  sum 
of  convective  and  conductance  matrix,  (S)  the 
source  vector,  and  { u)  the  time-dependent  temper¬ 
ature  values  at  nodal  points  of  the  solution 
domain.  A  Newton-Ralphson  iterative  procedure 
[6]  has  been  used  to  solve  the  equation  set  (4) 
for  the  {a(t,x))  values,  and  a  computer  program 
UDCFLOW  has  been  developed  to  facilitate  the 
calculations.  In  addition,  UDC  has  developed  an 
interface  between  the  PDA's  PATRAN  modeling  soft¬ 
ware  and  the  UDC's  UDCFLOW  finite  element 
analysis  computer  program. 

The  PATFJLN  program  provides  an  interactive, 
graphics-oriented  environment  for  creation  of  a 
geom.etric  and  finite-element  model  and  for  dis¬ 
play  of  analysis  results.  The  interface  allows 
PATRA-N  to  be  used  for  the  model  creation  and 
post -processing  steps  and  UDCFLOW  to  be  used  for 
the  analysis.  The  interface  consists  of  a  stand¬ 
alone  translator  program  PATUDC  that  will  read 
the  PATRAI'i  model  data  and  transform,  it  into  a 
UDCFLOW  input  file.  After  fluid/thermal  analysis 
by  UDCFLOW  is  complete,  the  results  can  be  trans¬ 
ferred  back  to  PATRAN  for  post-processing  dis¬ 
play. 

CLASSROOM  ADAPTATION  AND  WORKED  EXAMPLES 

For  classroom  instruction  of  the  thermal/fluid 
system  modeling  using  tne  PATRAN  and  UDCFLOW 
softwares,  students  m.ay  first  be  introduced  to 
the  use  oi  PATRAN  to  generate  geometric  and 
finite  element  models  for  a  variety  of 
f iuid/thermal  systems.  They  will  then  learn  to 
run  the  UDC  developed  PATRAN/UDCFLOW  interface 
computer  progra,mi  PATUDC  which  will  read  PATRAN 
model  data  and  transform,  them  into  UDCFLOW  input 
files.  Source  codes  for  a  version  of  the  UDCFLOW 
computer  program  will  then  be  given  to  the 
students.  After  the  students  have  familiarized 
themselves  with  the  UDCFLOW  program,  they  will  be 
assigned  to  solve  a  number  of  .problems  with  dif¬ 
ferent  degrees  of  com.plexity.  Etveral  examples 
of  tnese  problems  are  presented  below. 


However,  at  time  t=0,  the  electric  current  is 
shut  off,  thus  causing  a  transient  in  the  tem¬ 
perature,  while  the  slab  faces  continue  to  be 
held  at  100  degrees  C  by  a  coolant.  The  thermal 
diffusisity  of  the  material  is  constant  at  the 
value  of  2.84x10  'mVs.  We  would  like  to  find 
the  temperature  within  the  slab  as  a  function  of 
position  and  time  by  using  the  PATRAN/UDCFLOW 
softwares . 


A  model  for  this  problem  can  easily  be  created 
by  using  the  following  PATRAN  commands: 

GR,1,,0 
GR,2,TR,0.3,1 
LI, 1, ST, ,1,2 
GF,1L,G,33 
CF.IL.BAR 
PM, 1, TAN, . . . 

PH, 2, TAN, . . . 

PH, 3, TAN, . . . 

PM, 4, TAN, . . . 

PH ,  5 , TAN , . .  . 

PM, 6, TAN, . .  . 

DF,1L . 1 


The  PATRAN  output  file  PATRAN. OUT  can  then  be 
read  by  PATUDC  interface  program  to  produce  a 
UDCFLOW.DAT  file  which  can  be  read  by  the  UDCFLOW 
finite  element  analysis  program.  Output  file  from 
the  UDCFLOW  run,  UDCFLOW. RST,  can  then  be  used  by 
the  PATRAN  post-processor  P/PLOT  to  produce  th: 
X-Y  plot  ‘  temperature  versus  time.  Figure  1 
shows  temperature  distributions  of  the  slab  at 
time  equal  to  0,  150,  300,  450,  and  750  seconds, 

respectively. 

X  -  FCilTK-x 

3.  .350  .10  .:S  .iC  .86  .JC 


To  illustrate  the  procedure,  we  will  beqin 
with  a  sim.ple  1-D  transient  heat  conduction  prob¬ 
lem,  with  known,  boundary  temperatures.  This  will 
be  followed  by  more  difficult  transient  probiem.s 
which  begin  to  display  the  power  and  importance 
of  the  modeling,  analysis  and  display  softwares. 

Consider  an  infinite  slab  of  thickness  0.3 
meter  whicii  was  discussed  in  reference  (7). 
Initia’ly  because  of  an  electrical  current 
passing  the  slab,  the  slab  is  in  the  steady  state 
with  a  temtperature  distribution  in  degrees  C  i.s 
given  by  the  equation: 

T  -  100  ♦  400sin(-  x/0.3).  (5) 


Fig.  1  Simulated  Transient  Slab  Temperature 


For  problems  with  exitance  of  combined  convec¬ 
tive,  diffusion  and  internal  energy  source,  let 
us  consider,  for  example,  a  problem  which  has  the 
following  data: 

A  =  0.0013952  ft 
q  =  2,400  sin(  x/L) 
c  =  50  Ib/ft 
U  =  1.2  ft/s 
C  =  1.1  Btu/lb.F 

k  =  1.028  But/ft. F.s.  («■) 
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Initiall..  the  fluid  is  uniformly  at  200  degrees 
F.  At  t  greater  than  zero,  the  temperature  of  the 
fluid  is  held  at  20C  degrees  at  x  equal  to  zero, 
and  gradient  of  the  temperature  dT/d>:  is  main¬ 
tained  to  be  zero  at  x  equal  to  L.  Figure  2 
shows  resultant  temperature  distributions  of  the 
fluid  at  several  different  time  periods  of  t 
equal  to  1,  2,  3,  5,  and  9  seconds,  respectively. 


;  ‘  ?  c  ‘  e 


Fig.  2  Simulated  Transient  Fluid  Temperature 


CONCLUSIONS 


Computer-aided  modeling  of  thermal  and  fluid 
systems  has  gained  acceptance  both  in  scientific 
research  and  industrial  design.  Typical  CAD/CAE 
systems  for  f iuid/thermal  modeling  include  finite 
ele.ment  analysis  code  and  graphic  facilities  for 
pre-  and  post-processing.  To  instruct  the  stu¬ 
dents  early  in  their  career,  the  authors  have 
developed  at  UDC  a  finite  element  fluid/thermal 
analysis  computer  program  UDCFLOW  which  has  beer, 
integrated  with  the  PDA's  interactive  graphics- 
oriented  finite  element  modeling  pre-  and  post¬ 
processing  software  PATRAN.  It  resulted  in  an 
efficient  system  which  the  students  can  use  to 
gain  hands-on  experience  in  computer-aided 
modeling  of  the  fluid/thermal  systems  and  to 
solve  more  involved  problems  than  was  previously 
possible.  It  is  hoped  that  this  will  enhance  the 
students'  learning  experience  and  increase  their 
problem  solving  abilities. 
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